function [ifaceb,nodes] = pfbou2d(xnod,icone,ijob)
%
%    [ifaceb,nodes] = pfbou2d(xnod,icone,ijob)
%
%  If ijob = 1 --> calcula ademas del arreglo [elemento cara] 
%                  tambien el arreglo [node_1 node_2] de cada elemento 
%                  del contorno
%

tol = 1e-10;

icone0 = icone;

if nargin <3, 
   ijob=0;
   nodes=[];
end

% remuevo los elementos que tienen nodos duplicados (colapsados)
jcone = icone';
jcone = sort(jcone);
mask = sum(diff(jcone)==0);
eles_to_remain = find(mask==0);
eles_to_remove = find(mask>0);
icone(eles_to_remove,:)=[];

[nnod,ndm] = size(xnod);
[nele,nen] = size(icone);

%if ndm~=2, error(' ndm should be 2 ');end
if nen~=3 & nen~=4, error(' nen should be 3 or 4 ');end

if nen==4,
inoca = [1 2 ; 
         2 3 ;
         3 4 ;
         4 1 ];
elseif nen==3,
inoca = [1 2 ; 
         2 3 ;
         3 1 ];
else
error(' nen should be 3 or 4');
end

[nface,mface]=size(inoca);

inofa = (nnod+1)*ones(nele*nface,mface);
ifael  = [];

for k=1:nface,
 faces = (1:nface:nele*nface)'+k-1;
 ifael(1:nele,k) = faces;
 mm = mface;
 for m=1:mm,
  inofa(faces,m) = icone(:,inoca(k,m));
 end
end

% Remuevo caras duplicadas
% Ordeno nodos por cara por orden decreciente 
  [jnofa,is0] = sort(inofa');
   jnofa      = jnofa';

% ordeno por multiples llaves
  [knofa,is1] = mksort(jnofa);

% calculo operador inverso de ordenamiento
  ism1(is1,1) = (1:length(is1))'; 

% Reordeno el arreglo de caras por elemento
  for k=1:nface,
   ifael(:,k)= ism1(ifael(:,k));
  end

% Calculo aquellas caras duplicadas
  mask = (sum(abs((diff(knofa))'))'~=0);
  mask = [1;mask];

  jfaces = find(mask==0);
  jfaces = m2v([jfaces-1 , jfaces]);

  kfaces = (1:length(mask))';
  kfaces(jfaces) = zeros(length(jfaces),1);

clear jfaces mask knofa inofa jnofa ism1 is1

ifaceb = [];

  for k=1:nface,
   iaux = kfaces(ifael(:,k));
   iaux = (find(iaux~=0));
   ifaceb = [ifaceb ; [iaux ones(length(iaux),1)*k ] ];
  end
  
  if isempty(ifaceb)==1,      
     disp(' NO HAY ELEMENTOS EN EL CONTORNO ')
     nodes = [];
     return
  end

ifaceb(:,1)=eles_to_remain(ifaceb(:,1))';

if ijob==1,
faces = ifaceb(:,2);
eles  = ifaceb(:,1);

[numel,nen]=size(icone);
inod = faces+1;
mask = (inod>nen);
nn=find(mask==1);
inod(nn)=ones(length(nn),1);
inod = [faces,inod];

node_1=getcol(icone0(eles,:),inod(:,1));
node_2=getcol(icone0(eles,:),inod(:,2));

nodes =  [node_1,node_2];
end
